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ABSTRACT 


During this research period we have made significant progress in the four proposed 
areas: (1) Design of robust controllers via H“ optimization, (2) Design of robust controllers 
via mixed H 2 /H°° optimization, (3) M-A Structure and robust stability analysis for 
structured uncertainties, and (4) A study on controllability and observability of perturbed 
plant. 

It is well known now that the two-Riccati-equation solution to the H°° control 
problem can be used to characterize all possible stabilizing optimal or suboptimal H°° 
controllers if the optimal H°° norm or y, an upper bound of a suboptimal H°° norm, is 
given. In this research, we discovered some useful properties of these H°° Riccati 
solutions. Among them, the most prominent one is that the spectral radius of the product of 
these two Riccati solutions is a continuous, nonincreasing, convex function of y in the 
domain of interest. Based on these properties, quadratically convergent algorithms are 
developed to compute the optimal H°° norm. We also set up a detailed procedure for 
applying the H°° theory to robust control systems design. The relationship between the H“ 
norm and robustness issues has been carefully reviewed and the guidelines to formulate H°° 
optimization problems including the construction of a state-space realization of the 
generalized plant have been established. The controller formulas of Glover and Doyle are 
slightly modified and used to construct an optimal controller without any numerical 
difficulty. 

The desire to design controllers with H°° robustness but H 2 performance has 
recently resulted in mixed H 2 and H°° control problem formulation. The mixed H 2 /H°° 
problem have drawn attentions of many investigators. However, solution is only available 
for special cases of this problem. We formulated a relatively realistic control problem with 
H 2 performance index and H°° robustness constraint into a more general mixed H 2 /H°° 
problem. No optimal solution yet is available for this more general mixed H 2 /H°° 
problem. Although the optimal solution for this mixed H 2 /H°° control has not yet been 
found, we proposed a design approach which can be used through proper choice of the 
available design parameters to influence both robustness and performance. 

For a large class of linear time-invariant systems with real parametric perturbations, 
the coefficient vector of the characteristic polynomial is a multilinear function of the real 
parameter vector. Based on this multilinear mapping relationship together with the recent 
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developments for polytopic polynomials and parameter domain partition technique, we 
proposed an iterative algorithm for computing the real structured singular value. The 
algorithm requires neither frequency search nor Routh's array symbolic manipulations and 
allows the dependency among the elements of the parameter vector. Moreover, the number 
of the independent parameters in the parameter vector is not limited to three as is required 
by many existing structured singular value computation algorithms. 

For task 4, the work during this period concentrated on developing software for 
investigating the controllability robustness of linear time invariant systems. A preliminary 
software package (described in Appendix A) was developed for measuring the size, shape 
and location of the recovery region of initial states in finite time by bounded control. 
Although the various optimization algorithms (needed to automatically obtain the value of 
the indicators) work, they need to be refined in order to take care of such problems as ill- 
conditioned recovery regions. 
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PROGRESS REPORT 


ROBUST CONTROL OF SYSTEMS WITH REAL PARAMETER 
UNCERTAINTY AND UNMODELLED DYNAMICS 

1. INTRODUCTION 

This document is the second-period progress report on the NASA supported 
research, "Robust Control of Systems with Real Parameter Uncertainty and Unmodelled 
Dynamics", (No. NAG-1- 1102). We are happy to report that in this research period we 
have made significant progress in the following proposed research problems: (1) Design of 
robust controllers via H°° optimization, (2) Design of robust controllers via mixed H 2 /H°° 
optimization, (3) M-A Structure and robust stability analysis for structured uncertainties, 
and (4) A study on controllability and observability of perturbed plant 

Doyle, Glover, Khargonekar, and Francis (abbr.: DGKF) [1], and Glover and 
Doyle [2] presented a celebrating two-Riccati-equation type solution to a standard H“ 
control problem. The two-Riccati-equation method characterizes all possible stabilizing 
suboptimal H°° controllers whose order is not higher than that of the plant. The suboptimal 
H“ controller formulas in [1], [2] can be easily transformed into descriptor forms to 
construct optimal H°° controllers without numerical difficulties if the optimal H°° norm is 
given. The optimal H°° controllers, with very few exceptions, have direct feedthrough 
terms and therefore infinite bandwidth. Hence, control engineers may prefer strictly proper 
suboptimal H°° controllers to the optimal ones. However, knowing the optimal H“ norm is 
important in determining which suboptimal controller to be chosen in practical design. 

Recently, an efficient algorithm for computing the optimal H~ norm was proposed 
by Scherer [3]. Scherer considered the inverse (or pseudo inverse) of the DGKF H°° 
Riccati solutions, Xoo(y) and Y^fy), defined a new independent variable |i = y ", and 
showed that these inverses are concave functions of |i in matrix sense on their domains of 
definition. Based on this fact, a quadratically convergent Newton-like algorithm was 
proposed to compute the optimal H°° norm. 

Pandey et. al.'s hybrid gradient-bisection method [4] and Chang et. al.'s double 
secant and bisection method [5] were also proposed for the computation of the optimal H°° 
norm. The significance of the conjecture that p[X 00 (y)Y 0 .(y)], the spectral radius of 
X„(y)Y 00 (y), is a convex function of y 2 was mentioned in these two papers. Since there 
was no proof for this conjecture, bisection was used in these two algorithms as supplement 
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to guarantee convergence. 

In this research period, we discovered several important properties of the two 
DGKF H“ Riccati solutions, X 00 (y) and Y„(y). Among them, the most prominent one is 
that p[X 00 (y)Y 00 (y)] is a continuous, nonincreasing, convex function of y on ((3, °°), where 
P is the infimum of y such that the two DGKF H“ Riccati solutions, X^y) and Y„(y), exist 
and are positive semidefinite. Based on this property, a quadratically convergent algorithm 
can be easily developed to compute the optimal H°° norm, i.e., the such that 

p[X»o(Yoo)Yoo(Yoo)] = if a starting y is given inside the interval (P, i.e., y > P and 

p[X4y)Y 0 .(y)]>y 2 . 

A y inside the interval (P, y*,) most of the time can be easily found without the 
knowledge of p. However, the computation of p can be necessary if P itself is the 
optimum, which occurs when p[Xoo(p)Yco(P)] < p 2 . Newton-Raphson’s method can be 
employed to compute p and the convergence is quadratic. 

H°° control theory can handle the following two robustness issues: (1) Minimization 
of the maximum error energy for all command/disturbance inputs with bounded energy, 
and (2) Closed-loop stability under unstructured plant uncertainties with bounded H“ 
norm. Besides, H~ control theory is also indispensable in the structured singular value 
treatment of structured plant uncertainties [6,7], Detailed procedures to formulate these 
robust control problems into H°° optimization problems will be given in the report. 

For most practical control systems, a particular performance index is usually of 
considerable interest. At the same time there is also a desire to improve on the closed-loop 
system robustness. However, it is well known that there exists a tradeoff between these 
two design objective. Any improvement gained in one of the design goals is usually 
accompanied by a loss in the other. The Linear Quadratic Gaussian (LQG) Theory has 
been used successfully to design observer-based controllers with optimal performance for 
plants with fixed (or fixed power spectrum) exogenous signals. It is well known however, 
that the controllers derived using this approach possess poor robustness properties when 
compared to the guaranteed stability margins provided by full-state feedback control. 
Doyle and Stein [8, 9] devised a procedure, usually referred to in the literature as the Loop 
Transfer Recovery (LTR), to asymptotically recover the full-state feedback loop by tuning 
the LQG designed observer. The adjustment procedure is achieved by introducing a 
fictitious noise to the nominal plant input. A new LQG observer based controller is then 
derived. For the case of minimum-phase plants, the asymptotic recovery is achieved by 
letting the intensity of the added fictitious noise approach infinity. The procedure however. 
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has limitations when applied to nonminimum phase plants. 

In this research we concentrate on the recently developed H°° theory which evolved 
from the sensitivity minimization problem formulated by Zames [10]. It was shown there 
that many of the classical design objectives can be incorporated in the H°° design. For 
instance, the modeling of plant uncertainties can easily be formulated in terms of normed 
H°° plant neighborhoods. Another advantage of this relatively new theory can be revealed 
in the flexibility it offers to the designer in the formulation of the problem. For example, 
the well known design methodology of loop shaping can easily be achieved via the choice 
of frequency dependent weights on input as well as output signals. Although no direct 
relations between these two are known, with the new state-space solution to the general 
H°° problem [1, 2, 11] the design approach as will be shown later reduces to simply 
varying certain design parameters in order to achieve the design objectives. 

The motivation to design controllers for desired H 2 performance with H°° robustness 
constraints has recently resulted in mixed H 2 and H°° control problem formulation [11, 12, 
13, 14]. Special cases of this problem have been studied there and a solution in the form of 
coupled Riccati equations has been proposed. Although the optimal solution for this mixed 
H 2 /H°° control has not yet been found, we proposed a design approach which can be used 
through proper choice of the available design parameters to influence both robustness and 
performance. 

Stability robustness is an important issue in the analysis and design of control 
systems. Currently, there are two major approaches to stability robustness analysis. One is 
the structured-singular-value (SSV) [6,7] or the multivariable-stability-margin (MSM) 
[15,16] approach and the other is the perturbed-characteristic-polynomials approach [17- 
22]. Several significant progresses have been made in both approaches. In this research, an 
iterative algorithm of computing the real SSV and the real MSM is developed based on the 
existing results in both approaches. 

For a large class of linear time-invariant systems with real parametric perturbations, 
the coefficient vector of the characteristic polynomial is a multilinear function of the real 
parameter vector. Based on this multilinear mapping together with the recent results by De 
Gaston and Safonov [15], Sideris and Pena [16], Bartlett, Hollot, and Lin [18], and 
Bouguerra, Chang, Yeh, and Banda [23], an algorithm for computing the real structured 
singular value is proposed. The algorithm requires neither frequency search nor Routh's 
array symbolic manipulations and allows the dependency among the elements of the 
parameter vector. Moreover, the number of the independent parameters in the parameter 
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vector is not limited to three, as is required by many existing structured singular value 
computation algorithms. 

The literature on controllability gives various degree of controllability (DOC) 
measures, each based on different point of view. That is either, in terms of the 
controllability grammian or the recovery region of initial states in finite time by bounded 
control, or the variation of the system and control matrices. The work during this period 
concentrated on: (i) trying to develop a relation between the various DOC measures 
encountered in the literature, and (ii) developing software for investigating the 
controllability robustness of linear time invariant systems. A preliminary software package 
(described in Appendix A) was developed for measuring the size, shape and location of the 
recovery region. Although the various optimization algorithms (needed to automatically 
obtain the value of the indicators) work, they need to be refined in order to take care of 
such problems as ill-conditioned recovery regions. 

In section 2 of this report, we will show the overall progress in this research 
period. The work performed and the status of the proposed tasks are summarized in section 
3. Section 4 is the conclusion. The work for future research will also be briefly described 
in section 4. 


2. OVERALL PROGRESS 


Consider the system 


z(s) 

y(s) 


G n (s) G 12 (s) 
G 21 (s) G 22 (s) 


v(s) 

u(s) 


:= G(s) 


v(s) 

u(s) 


(2- la) 


u(s) = K(s) y(s) (2- lb) 

where G n (s) e R(s) PlXmi , G 12 (s) e R(s) PlXm2 , G 21 (s) e R(s) P2Xm ', and G 22 (s) e 
R(s) P2Xm2 . IR(s) pxm is the set of pxq proper rational matrices with real coefficients. Recall 

that the standard H°° optimization problem is the problem of finding a proper controller 
K(s) such that the closed-loop system is internally stable and IIT^II^ is minimized where 

T zv (s) is the transfer function of the closed-loop system from v to z. 


The realization of the generalized plant G(s) is given by 
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G(s) = 


with the following assumptions: 

(i) (£j. A) has no unobservable modes on jto-axis and (A, B t ) has no uncontrollable 
modes on j(i)-axis, where 


A — A-B 2 D n C u 

(2- 3a) 

^i= 

(2-3b) 

A = A-B^^Cj, 

(2-3c) 

B i = B r B i D 2i D 2i- 

(2-3d) 

0, D 12 = j , D 21 = [0 I], d 22 = 0. 



(iii) (A, B 2 ) is stabilizable and (C 2 , A) is detectable. 

The two Riccati equations involved in the H“ solution can be expressed as: 

A T + X„ A + X M ( y 2 BjBj 1 - B 2 B^ ) = 0 (2-4a) 

and 

AY 00 + Y„A t +Y 00 (y 2 c]'Ci-C^C 2 ) Y^ + BiBj = 0, (2-4b) 

where A, A, Bi and are defined on (2-3). 

The following theorem characterizes suboptimal stabilizing controllers such that IIT^IL < y. 
DGKF Theorem : [1] 

There exists a stabilizing controller such that IIT^IL < y if and only if the following 
three conditions hold. 

(i) There exists a positive semidefinite stabilizing solution X m (y) to (2 -4a). (2-5a) 

(ii) There exists a positive semidefinite stabilizing solution Y„(y) to (2-4b). (2-5b) 

(iii) p[XJy)Y„(Y)]<Y 2 - (2* 5c) 
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Moreover, when these conditions hold, one such controller is 


K sub( S > = 

r\ 

Bk_ i 

_ c k 

0 

with 



B k = -EL 2 , 


C k =F 2 


( 2 - 6 ) 


Ak 


A + 


B , B 2 


+ EL 2 (C2+D2iFj), 


where 

E = (I- y 2 YJ„ ) 1 , F, = -Bfx.., F 2 = -B T 2 X„- D^Cj, 

Lj =Y 2 Y 00 C{, L2 = -Y^Cj- BjDjj. 


The above theorem shows an easy state-space approach to construct a stabilizing 
suboptimal controller such that IIT ZV IL < y. The theorem can also be used to compute the 
optimal H“ norm and construct an optimal H” controller. The optimal IIT^IL is the infimum 
of y such that the above three conditions hold. With very few exceptions, the optimum 
occurs when ptX^^Y^y)] = Y 2 which will render A k and E undefined. This numerical 
difficulty can be easily resolved by transforming (2-6) into a descriptor form [1][25][26]. 


2.1 PROPERTIES OF H~ RICCATI SOLUTIONS 

First, we assume (Cj, -A) is detectable. This assumption will be removed later in this 
section. Suppose we have Riccati equation 

&A t + A& + (Y 2 BiB{- B 2 B 2 ) + =0, (2-1) 

then solution &(y) has following properties. 

Theorem 2.1: &(y) is a well-defined function on (a*, +<»), where 

a x := inf {y : Ye R+ and ft(y) exists} (2-8) 

Moreover, &(y) is analytic with 

hi) := ^(Y)) > 0 (2-9) 
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( 2 - 10 ) 


&(Y) := -S$(Y)) * 0 

dr 

on (a x , +oo). 

Furthermore, the following theorem gives the eigen properties of &(y). 

Theorem 2.2: On (a x , +«>), 

a) all eigenvalues of X(y) are smooth, nondecreasing functions of y ; 

b) the minimal eigenvalue of &(y) is a nondecreasing, concave function of y ; 

c) ft(y) is invertible almost everywhere. 

This theorem implies that if we define 

(3 X := inf {y: ye 1R+ and ft(y) is positive semidefinite}, (2-11) 

then either A. min [&(P x )] = 0 or P x = a x , which in turn implies &(y) is positive definite 
everywhere on (|3 X , +<»). 

With these properties of ft(y) in mind, one can easily find the corresponding 
properties of X^y) by comparing the following two equations. 

A t X 00 +X 00 A+XJy 2 B 1 b}-B 2 b5x oo +c]'Ci = 0 (2-12) 

&A t + A& + (y 2 BiB{- B 2 B 2 ) + = 0. (2-13) 

In the beginning of the section, we assumed that (Cj, -A) is detectable. In this case, it is 
easy to see that X oo (y)= & _1 (y) almost everywhere on (a x , +««) and X oo (y)>0 on (|3 X , +<*>). 
This assumption will be removed in the following. 

In the case that (Cj, -A) is not detectable, one can always find an orthogonal matrix 

U = [Ui U 2 ] (2-14) 

such that 

UjAU, 0 

UjAU, U^AII, 

1 1 



u t au = 


(2-15) 



r t t ”i 
U>1 u 1 B 2 

U t B = [u T B. U T Bj = _ T (2-16) 

' U‘B,_ 

and 

CjU-lCiU! 0] (2-17) 

with (CiUi, -U^AUj) detectable [27]. Furthermore, the solution to (2-12) can be 
expressed as 

" X(y) 0] 

X~(Y) = U[ 0 0 JU T (2-18) 

with X(y)= & _1 (y) almost everywhere on (a x , +»), where X(y) and 5((y) are the stabilizing 
solutions to (2-12) and (2-13) respectively with (A, Bi, B 2 , Cj) replaced by (uJaUj, 
T T 

UjBj, U t B 2 , QUO. Therefore, no matter whether (C p -A) is detectable or not, it is 

always true that X m (y) exists almost everywhere on (a x , +°°) and X oo (y) > 0 on (P x , +°°). 

Moreover, we have the following theorem which gives the properties of the first and 
second derivatives of X m (y). 


Theorem 2.3: On (P x , + 00 ), 


a) 

X_(Y) = ^X_(Y)> £ 0, 

(2-19) 

b) 

X_(Y) ^X-(Y» a 0. 

(2-20) 


Again, next theorem gives the eigen properties for Xjty). 

Theorem 2.4: On (P x , +°°), 

a) all eigenvalues of X oo ( y) are smooth, nonincreasing functions of y ; 

b) the maximal eigenvalue of X oo (y) is a nonincreasing, convex function of y. 

Similarly, we can obtain the same results for Yjy), the stabilizing solution to 

Y„A t + AY^ + Y oo (y’ 2 c}c 1 - c£ 2 )Y 00 +B 1 b'[ = 0. (2-21) 
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If (-A, BO is stabilizable, then Y Jy)= ^ l (y) almost everywhere on (oty, +°°) and 
YJy)> 0 on ((3 y , -h»), where ^(y) is the stabilizing solution to 

A T t + tA + (y' : 2 cjc r C^+tBiB^ = 0 (2-22) 

and 

oty := inf {y : ye 1R + and Y^y) exists} (2-23) 

P y := inf {y : ye 1R+ and Y m (y) is positive semidefinite } . (2-24) 

If (-A, Bi) is not stabilizable, then an orthogonal matrix V= [Vi V 2 ] can be found such 
that 


V 1 AV 1 V 1 AV 2 


v t av= 


0 


V^AV,, 



T 

V B, 


' C i V “ 


C V 

H V 1 

c V 

S 2 

, v T B,= 

1 1 

_ 0 . 

and CV = 

_ C 2 V _ 

— 

C V 

L 2 v 1 

C V 
v_ 2 v 2 J 


T T 

with (-Vj AV p VjBi) stabilizable. Furthermore, the stabilizing solution to (2.21) can be 
expressed as 


Y«(y) = v 


Y(y) 0 

0 0 


(2-25) 


with Y(y)= ^ _1 (y) almost everywhere on (a y , <») and Y(y)> 0 on (P y , +<»). Again, Y(y) 
and ^(y) are the stabilizing solutions to (2-21) and (2-22) respectively with (A, Bi, Ci, 
T T 

C 2 ) replaced by (VjAVj, VjBi, CjVi, C 2 V 1 ). Note that all results presented above for 

Xooly), X(y) and &(y) have their counterparts for Y«(y), Y(y) and V(y) respectively. 

Before moving to the next theorem to investigate the properties of X^Y*,, we define a 
and P as follows: 

a := max{oc x , oty } ( 2 - 26 ) 

P := max(P x , p y }. (2-27) 

With these definitions, we can see that X»Y« exists on (a, +°°) almost everywhere. 
Moreover, X^Y.* has no negative eigenvalues on (P, +°°), since both X„ and Y« are 
positive semidefinite on (P, +■»). Now, we are in the position to state our main result. 


Theorem 2.5: On (P, +°°), 

a) all eigenvalues of X^Y^ are smooth, nonincreasing functions of y ; 

b) pCXooYoo) is a nonincreasing, convex function of y. 


1 3 



Furthermore, if we define p(y):= p(XcoYoo) and p(y) 


d(p(y)) 

dy 


on (p, +»), the slope 


of p(y) can be expressed as 


p(y) = 


v T (X 00 Y 00 +X 00 Y 00 )w 


T 

V w 


(2-28) 


where v and w are the left and the right eigenvectors of X^Y,,, respectively corresponding 


to its maximal eigenvalue. X,,, and Y^ can be obtained by solving the following Lyapunov 
equations: 


A T X m +xJ. -2y 3 X oo B 1 B[X oo = 0 


(2-29) 


and AY^+'La 1 ’ -2y 3 YX]'C 1 Y 00 = 0, (2-30) 

with A = A + (y 2 BiB|-B 2 B 2 )X oo and A = A + Y^cfo-C^C;,). 


2.2 ALGORITHMS TO COMPUTE THE OPTIMAL H~ NORM 

According to DGKF theorem, we can see that finding the optimal H M norm, 
denoted as y„, is equivalent to finding the infimum y such that all three conditions in (2-5) 
hold. From the previous subsection, it is obvious that y„e [P, +<»). It is possible for P to 
be y„, especially when P and a are identical, however, with very few exceptions, y^e (P, 
+<*>), which implies that y^ is the solution to p(y) = y 2 . The relations between a, P and y M 
are shown in the figure below. 

Both X~ and Y°» exist 

Both X~ > 0 and Y~ > 0 


All three conditions hold 



0 a P 

Fig. 2.1 


The figure implies that the problem of finding the optimal y„ is actually that of either 
searching for the intersection point of p(y) with y 2 inside (P, +«>) or computing the 
boundary point p. Since p(y) is a convex function on (P, +■»), then gradient searching 
method can be employed. 

Assume that we have a starting point y n in the interval (P, y„), then the optimal y 
can be obtained easily as follows. Refer to Fig. 2.2, draw the tangent line with slope p(y n ) 
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passing through the point (y n , p(y n )). The abscissa, y n+1 , of the intersection of the tangent 
line and the curve y = y 2 , always lies between y n and y*. The search process is repeated 
until the gap y*, - y n+1 is small enough. 



Furthermore, we will see that the convergence rate is quadratic. Define e n = y„ - y n and 
e n+ i = y„ - y n+1 . It is straightforward to show that 


e n+l ~ 


P(YJ 2 

2 I p(YJ - 2y M I 6,1 


(2-31) 


which implies quadratic convergence. For convenience, the algorithm described will be 
referred to as Q-step, since the quadratic convergence is guaranteed, provided there is a 
starting point y n e (P, yj to start with. 

Based on the discussion above, following procedure is given to compute the 
optimal norm. 

Step 1 . Initial point 

Refer to Fig. 2.1, choose a relatively large initial y x such that Yi€[| 3, +°°). If 
Yie(p, yj, then it can be used as the starting point for y^. Go to Q-step. If (y«,* 
+»), which implies X^Xyj) > 0 and > 0, then go to step 2. 

Step 2. Starting point for y„ 

Evaluate p(y) at y L , we have the point (yj, p(Yi)). Refer to Fig. 2.3, draw a line 
passing through the point (yj , p (Yi)) with slope p(Yi). The abscissa, Y 2 , of the 
intersection of the straight line and the curve y = y 2 , is always less than y„. If y 2 e (P, 
y^), then we are ready to go to Q-step with y 2 as the starting point y n in Fig. 2.2. 
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Remark: Most of the time, the method described in the previous paragraph gives the 

starting point y n inside the interval ((3, y^,). However, this method may fail, see Fig. 2.4 
and Fig. 2.5. For the case described by Fig. 2.4, a smaller Yj could do the job. However, 
for the case of Fig. 2.5, y 2 is always less than a, and therefore less than (3, no matter how 
small Yj is. Since it is difficult to tell which case we are facing and there is no efficient 
guideline to reduce y t , we suggest to continue if two or three rials of Yi does not give a 
starting point. 

If y 2 € (ot, (3), go to step 3 to compute (3. If y 2 < <*, then go to step 4 to compute a. 



Step 3. Computing (3 

For the computation of [3, one also needs a starting point to use Newton method. It 
is similar to the gradient method described earlier. See Appendix for details. Note that 
the Y 2 obtained in step 2 can be used as a starting point for |3. If |3 is optimal y (i.e., all 
the three conditions in (2-5) are satisfied), see Fig. 2.1, then the algorithm will 
terminate here. Otherwise (3+e can be served as the starting point y n for y^ in Fig. 2.2, 
where e is a very small positive real number. Go to Q-step. 

Step 4. Computing a 

If either X 00 (y 2 ) or Y oe (y 2 ) does not exist, which implies y 2 < a, we will compute a 
first. Appendix will give the algorithm for the computation of a. If a is optimal y (i.e.. 
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all the three conditions in (2-5) are satisfied), see Fig. 2.1, then the algorithm will 
terminate here. Otherwise a can be used either as a starting point for (3, when a*(3; or 
a starting point for y„, when a = (3. Go to step 3 and Q-step respectively. 

Q-step. Computing y^ 

This step was described in Fig. 2.2. We can see from the earlier discussion, once 
we have a starting point for y„, the quadratic convergence is guaranteed. Algorithm 
terminates. 

Example 2.1: The following is a simple H“ optimization problem which is used to 

illustrate the proposed algorithm of computing the optimal H“ norm. A realization of the 
generalized plant G(s) is given by 



' A 

B, 

B~ ‘ 


'-i 

0 

1 

0 

O' 



l 

2 


0 

2 

0 

0 

1 

G(s) = 


D n 

D 12 

= 

i 

1 

0 

0 

0 



D 

D 


0 

0 

0 

0 

1 


c 




wj-uwitu 




2 

21 

22_ 


1 

1 

0 

1 

0 


Starting from y 0 = 100, we have p(y 0 ) = p[X«,(Yo)Y«(yo)] = 2.16e01. The slope of 
p(y) at this point (y 0 , p(Yo)) is p(Yo) = -3.49e-05. The tangent line at this point will 
intersect with the curve y =y 2 at y x = 4.65. Again, evaluate p(yj) = 2.24e01 and compute 

the slope p(y) at (y l5 p(Yi)). Since p(y t ) > yj\ y l is inside the interval (p.y^) and therefore 

from now on the convergence is guaranteed. The process is repeated until the gap between 

Vp(Yo) and Y n is small enough. The following data show that only four iterations are 
needed to reach the optimum, y op = 4.734160476390407, with accuracy better than 10 14 . 


Iter 

Yn 

J P(Yn) * Yn 

P(Yn> 

ru 

100 

9.535e+01 

-3.487e-05 

i 

4.64799876 1 538930e+00 

8.403e-01 


2 

4.734064423866624e+00 

9.440e-04 

-3.595e-01 

3 

4.734 160476276923e+00 

1.115e-09 

-3.595e-01 

4 

4.734 160476390407e+00 

7.105e-15 



By the formulas in (2-6) and the descriptor-form technique described in [26], we are 
able to construct an optimal controller as follows: 
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Kopt(s) 


-.87542 

-0.13925 

4.42042 

-4.73416 


with the H°° norm of the closed-loop system equals y op . Note that the optimal controller has 
a direct feedthrough term and thus has infinite bandwidth. If we choose y = 4.8 which is 
about 1.4% higher than y op , we have a suboptimal controller 


-8.67072e-01 

1.32928e-01 

-1.38959e-01 

-1.38320e+01 

-1.52323e+02 

4.73733e+00 

_ -9.30025e+00 

-1.49792e+02 

0 


which has a reasonable bandwidth and the closed-loop H°°norm, HT zv ll oo < 4.8 which is 
only 1.4% away from the optimal H°°norm. 

2.3 Formulation of H°° Optimization Problems 

Many control problems can be formulated as the standard H“ optimization problem. 
For the purpose of demonstration, two examples are given in the following. The first is a 
mixed-sensitivity optimization problem to be formulated as a two-block H“ optimization 
problem; the second is a disturbance reduction problem with measurement noise which 
turns out to be a four-block problem. 

A. Mixed-Sensitivity Optimization Problem 

Consider the following system: 

y(s) = P(s)u(s) + v(s) (2-32a) 

u(s) = K(s)y(s) (2-32b) 

where v(s) is disturbance, y(s) is output and K(s) is controller to be designed. It is well 
known that a smaller 11(1- PK) *11^ means a better disturbance attenuation, whereas a smaller 
NPKfl-PK)' 1 II ^ implies a better robust stability. Unfortunately, the H~ norms of (I-PK) 1 

and PK(I-PK) 1 may not be made small at the same time. If we make one of them smaller 
then the other will become larger. To have a trade-off between these two quantities, 
Kwakemaak [24] formulated the mixed-sensitivity problem as the problem of finding a 
controller K(s) which stabilizes the closed-loop system and minimizes NOIL where is 
given by 



d> = 


(2-33) 


Wj(I-PK) 1 
W 2 PK(I-PK) 

W j and W 2 are weighting matrices chosen by the designer according to the concrete 

situation. In other words, they depend on the characters of the disturbances and system 
uncertainties. Usually, the disturbances occur most likely at low frequency, therefore 
Wj(s) is chosen to be a low-pass filter to emphasize the error energy at low frequency. The 
plant uncertainty is also frequency-dependent; the higher the frequency is, the larger the 
uncertainties become. Hence, W 2 (s) is usually chosen to be a improper transfer function 
(but W 2 P(s) has to be a proper transfer function), which is analytic in closed right half 
plane. In the following, we assume that W^s) is strictly proper, W 2 (s) is a polynomial 
such that W 2 P(s) remains proper and both of them are analytic in closed right half plane. 

The problem of finding a K(s) which stabilizes the closed-loop system and 
minimizes llthll^ can be rearranged into the standard H°° optimization problem. Consider the 

following system: 


11 

N- -T 

W l 

0 

W t P 

W 2 P 

[:] 

yj 

. I 

p . 



u = Ky (2-34b) 

It is easy to show that the matrix d> defined by (2-33) is just the transfer function from v to 
[Zj T z 2 t ] t of the closed-loop system (2-34). Comparing (2-34a) with (2-la), we can see 

that 


G 


11 




W t P 

w 2 p 




= P. 


If P, W 2 P, and Wj have state-space realizations as follows 



V 

. c p 

d p_ 


W 2 P = 


_a p 

V 

_Cw2 

D w2_ 


W. 


Awl 

B W 1 

_Qwl 

D w i_ 


(2-35) 


(2-36) 


Then the generalized plant G(s) has a state space realization as shown in (2-2) with 
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0 

, B, = 

’ 0 ' 

, B. = 

" B P _ 

_ B W, C P 

A , 
wlj 

’ l 

B , 

L wl J 

’ 2 

_ B wl D P_ 


X, c p 

C wl" 

, D, , = 

D i" 

wl 

, = 

D ,D d " 

wl P 

1 

n 

to 

0 

’ n 

. 0 . 

’ 12 

1 

$ 

Q 

1 


c, = [ C P °] , D 21 = I , D 22 = Dp (2-37) 

Note that because W 2 is a polynomial, the A-matrix of W 2 P is same as that of P. 

B. Disturbance Reduction Problem 



Fig.2.6 A disturbance attenuation problem 

Consider the feedback system shown in Fig.2.6. P(s) is a given plant, W^s), 
i=l,2,3,4 are weighting matrices, and K(s) is the controller to be designed. The 
disturbance and noise are the outputs of W 3 and W 4 driven by d and n respectively. z l is 
the weighted error response; z 2 is the weighted control input. Let z T = [Zj T z 2 T ] T , v T = [d T 

n T ] T and assume that v is unknown but with its energy bounded by unity. The objective is 
to find a controller K(s) which stabilizes the closed-loop system and minimizes the worst 
llzll 2 , i.e., minimizes the H“ norm of T^, the closed-loop transfer function from v to z. T zv 
is given by 


T 


ZV 


WjPfl-KPj'Vj W 1 PK(I-PK)' 1 W 4 
W 2 KP(I-KP) 1 W 3 W 2 K(I-PK)‘w 4 


(2-38) 


Note that W 1 PK(I-PK) 1 W 4 and W 2 KP(I-KP) 1 W 3 ) are the output and input 
complementary sensitivity functions. Their H~ norms indicate the stability robustness of 
the closed-loop system for the multiplicative plant uncertainty introduced at the output and 
input respectively. W 2 K(I-PK) 1 W 4 is the control complementary sensitivity function 
whose H~ norm indicates the stability robustness of the closed-loop system for additive 
plant uncertainty. Hence, reducing IIT^II*, will also improve the robust stability of the 

closed-loop system. 
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It is easy to verify that the generalized plant of the system can be expressed as: 


v 


•w lP w 3 

0 

w,p- 

Z 2 

= 

0 

0 

W 2 

_y _ 


- PW 3 

W 4 

P . 


d 

n 


(2-39) 


That is. 


’ll 



’ WjPW 3 0 ' 


~ WjP" 


. 0 0 _ 

• G 12 = 

w 

L 2 J 


G 21 = [ pw 3 w 4 ] , o M = p 


22 


(2-40) 


If P, Wp i=l,2,3,4 have state-space realizations as follows: 


P = 


_A P 

B P _ 

_ C P 

D P. 


W = 


i \vi 

B wi 

c 

^Wl 

D W1 . 


i= 1,2, 3,4 


(2-41) 


Then the generalized plant G(s) has a state space realization as shown in (2-2) with 



Ap 0 0 B p C wJ 0 1 


' B P D W, 0 ' 


r b p 1 

A = 

B w,Cp A wl 0 B w| D p C wJ 0 

Bi = 

B M D P D »3 0 

b 2 = 

B wl D P 


O 

o 

«? 

o 

o 


o 

o 


B w2 


0 0 o a„ 3 0 


B , 0 

w3 


0 


— 1 

<* 

o 

o 

o 

o 


. 0 b m _ 

y 

. 0 . 


C, = 


D w,Cp C W1 o d w1 d„c w3 0 
0 0 C , 0 0 

w2 


D 


11 


D wl D P D w3 

o" 

D,. = 

~D ,D“ 

wl P 

. 0 

0. 

12 

* 

D , 

L w2 J 


C 2 = [^TP ® ^ ^P G w3 ^w4], D 2 1 = P"V^w3 D w 4 ], D 22 = D p . 


(2-42) 


Above { A,B,C,D } is the state-space representation for the generalized plant G(s). 


2.4 A Design Approach to Achieve H 2 /H°° Objectives 

The mixed H 2 and H°° control problem in its most general form has not been solved 
as of yet. However, taking advantage of the recent advances in H°° theory, we believe 
that it is possible to develop a simple design procedure that will address the performance- 



robustness problem. In the following, we consider the disturbance attenuation problem 
with control weighting as shown in Figure 2.6. The nominal plant is denoted by <P(s) and 
the energy bounded exogenous signals consist of plant input disturbances d and 
measurement noise ru The weighted error response and the weighted control input consist 
of z\ and respectively. 

Following the approach presented by Doyle et. al. in [1,2], this problem can easily 
be transformed to the general H 2 /H°° configuration. The task is then to design a stabilizing 
H°° controller K(s) such that the H 00 norm of the transfer function matrix from the 
exogenous input vector w T =[ d T n T ] T to the output vector z T =[ zi T z 2 T ] T is minimized. 
From Figure 2.6, this transfer function matrix is given by: 

Wj P(I-KP ) 1 W 3 WjPK(I-PK ) 1 W 4 
W 2 K P ( I - K P y 1 W 3 W 2 K ( I - P K ) 1 W 4 

where Wj (s), i = 1, ... ,4 are weights to be chosen by the designer appropriate to the plant 
and design objectives being considered. Usually, W3 is a low frequency filter indicating 
that the disturbances introduced at the plant input are low frequency signals. The rest are 
usually assumed to be high frequency filters to include measurement noise, unmodeled 
dynamics, and any other plant uncertainties that may occur at high frequencies. Notice that 
W 1 can be nonproper, hence providing the flexibility to consider frequency dependent 
output responses. In this note however, we will only consider weights that lead to 
controllers having the same order as that of the nominal plant 

The above disturbance rejection problem turns out to be a four-block problem and 
its solution is described in [2]. It has been realized however, that suboptimal controller, 
i.e. controllers satisfying II T zw lloo < y for some positive y greater than the optimal H°°- 
norm of the closed-loop transfer matrix, are much easier to characterize than optimal ones. 
It also turns out that the variable y can be used as a design parameter. This can easily be 
verified when the general mixed H-/H°° control problem is restricted to the case where the 
H 2 closed-loop transfer matrix is the same as that of the H°° problem. In this case, the 
suboptimal H°° controller approaches the H 2 /LQG controller as y approaches <». 

From the input/output relation described above, it can be seen that it is possible to 
address the performance/robustness problem via the design weights. For instance, the 
complementary sensitivity function which is a measure of the closed-loop robustness is 
represented by the (1,2) block of the closed-loop transfer matrix above with the weights 
W 1 and W 4 included. Thus by appropriately choosing these weights one can hope to gain 
on performance without giving up too much robustness. As we will see later in the 
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(2-43) 




example, this procedure when used effectively is comparable to the LTR technique as far as 
H 2 performance and robustness are concerned. 

Example 2.2 

Consider the following typical LQG problem: 



y = [ 2 1 ] x + n 


with E(d) - E(n) - 0; E[d(t)d(x)] = E[n(t)n(x)]= 8 (t-x). The following performance index 
is of particular interest 

J = I ( x T Q T Q x + u 2 ) dt where Q = 4 V5 [V35 1 ] 

Jo 

In the H°° formulation, we construct the generalized plant 



by considering the same nominal plant as in the LQG problem with the following weights: 


^ ttl 8 + s ^+ 2 + ~ anc * ^ 2 -3,4 ( s ) = a 2 , 3,4 where the ai’s are real. The optimal 


LQG performance index turns out to be equal to 493. Introducing a fictitious noise at the 
plant input as described by the LTR technique improves the stability margins of the closed- 
loop system, but at the expense of performance. The solid-line curve 3 in Figure 2.7 
represents the singular values of the complementary sensitivity function with the LQG 
controller implemented while curve 4 corresponds to the LTR design. The improvement in 
robustness using the LTR technique results in a loss in performance as summarized in the 
table below. The dashed line Curve 1 in Figure 2.7 corresponds to the implementation of 
the (a-design) H°° controller with the constants ai = <X 2 = a .3 = 1 and a 4 = 0 . 01 . The 
reason why 014 was varied was because it is directly related to the complementary 
sensitivity function as indicated by the input-output relations above. The suboptimal H°° 
controller was obtained for a value of ythat is 15% higher than the optimal H°° norm. 
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Figure 2.7 Complementary sensitivity function (db) 


Curve 2 corresponds to the implementation of a H°° controller (p-design) with the 
following data: 


Bi = pi 


35 

.-61 


0 

0. 


Ct = p 2 


20 VT 

. 0 


4V5" 

0 


Di2 = P 3 


0 

. 1 . 


D21 = p4[ 0 l] 


Notice that the case when all the constants Pi's equal to 1 corresponds to the special case 
of the mixed H 2 /H°° control where both the H 2 and H°° closed-loop transfer matrices are 
the same. Curve 2 which was obtained by setting P 2 and P 4 to be equal to 0.1. Notice that 
at high frequencies curve 2 is below curve 4 (LQG/LTR) indicating a better robustness 
measure. The resulting performance index in this case is the same as the LTR case, i.e. 
equal to 587. One interesting point can be seen in here by examining the open-loop 
Nyquist plots of the LQG/LTR controller and the H°° controller shown in Figure 2.8. The 
dashed-line curve which corresponds to the H°° controller implementation illustrates a 
phase margin that is almost 10 degrees better than that of the LTR design. There is 
however a loss of 1 db in gain margin. It is not clear exactly why this happens but we 
intend to look into the problem in further research. 


2 4 




1 



-5 -3 -1 

Figure 2.8. Nyquist Plots 


2.5 Computation of the Real Structured Singular Value 

For a large class of linear time-invariant systems with real parametric perturbations, 
the coefficient vector of the characteristic polynomial is a multilinear function of the real 
parameter vector. Based on this multilinear mapping together with the recent results by De 
Gaston and Safonov [15], Sideris and Pena [16], Bartlett, Hollot, and Lin [18], and 
Bouguerra, Chang, Yeh, and Banda [23], an algorithm for computing the real structured 
singular value is proposed. The algorithm requires neither frequency search nor Routh's 
array symbolic manipulations and allows the dependency among the elements of the 
parameter vector. Moreover, the number of the independent parameters in the parameter 
vector is not limited to three, as is required by many existing structured singular value 
computation algorithms. 



Fig.2.9 Standard structure for a perturbed closed-loop system. 

All the plant uncertainties, structured or unstructured, unmodeled dynamics or 
parametric perturbations, can be described by the block diagram shown in Fig.2.9. In 
Fig.2.9, A(s) = block diag { Aj(s), A 2 (s), ..., A m (s) } and M(s) is the nominal system 

which includes the nominal plant and the stabilizing controller. The structured singular 
value (SSV) is defined based on the above perturbation structure. The SSV is 
nonconservative scalar stability- margin measures for multivariable systems. 

Algorithms [6,7] to compute the SSV are available only for those cases where the 
number of perturbation blocks are less than or equal to three. The computational problem 
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for the cases with more than three perturbation blocks is still an unsolved problem. 

One important special case of plant uncertainties is the real parametric perturbation. 
In this case the perturbation matrix A(s) is a real diagonal matrix. The SSV defined for this 

case is called the real SSV. An iterative algorithm of computing the real SSV for real 
diagonal A was developed by De Gaston and Safonov [15] and generalized by Pena and 
Sideris [16]. There is no limitation on the number of perturbation parameters. However, 
this iterative algorithm is complicated since an extensive frequency search is required. 

In the following we assume that the perturbation matrix A in Fig.2.9 is real 
diagonal, i.e., A = diag { 8j, § 2 , ..., 8 m } and the nominal system M(s) is a rational matrix 

with real coefficients. If the parameters vary independently and - 1 < 8 . < 1 , i = 1 , 2 ,. ..,m, 

the parameter perturbation domain JS) can be described as a hyper-cube <0 with 2 m vertices 
(±1, ..., ±1) in the m-dimensional real space. In general, the perturbation matrix A can be 
written as A = diag {5,I ml , S 2 I m2 , ..., S r I mr } where I mi is the identity matrix of order mi 

and ml+m2+ ...+mr = m. That is, 8 =8,= =8 =8,, 8 , =8 , = - =8 =8 

etc. In this case, the parameter perturbation domain JO is an r-dimensional hyperplane 
inside the m-dimensional hypercube j$). The system is said to be robustly stable in 49 if 
and only if it is stable for every parameter vector 8 = [ 8j 8 2 ... 8 m ] T in 49 . Throughout 

the report, we may use "49 is stable" to replace "the system is robustly stable in 49" 
whenever there is no confusion. 

The real multivariable stability margin (real MSM) k M is defined as the largest real 

constant k such that the closed-loop system remains robustly stable in k49 where k<© is the 
enlarged (or shrunk) parameter perturbation domain of 49, i.e., 

k<£> := { 8: 8= [8 1 ,..,8 1 , 8 2 ..,8 2 , ..., 8 r ,..,8 r ] e R m 

and I 8j I < k, i=l,2,...,r } (2-44) 

The enlarged (or shrunk) hypercube of j$>, kJ0 is 

kj0 := { 8 : 8 € R m and I 8 ; I < k, i=l,2,...,m }. (2-45) 

Recall that the real structured singular value (real SSV) ji is defined as 

[i := [ min { k I det [I+M(j(o)A] = 0 for some 0) and A e X (k) } ] _1 (2-46) 

where 

X (k) = { A I diag {8jl ml , 8 2 I m2 , ..., 8 f I mr ] with I 8 . 1 < k, for all i } (2-47) 

That is, the real SSV p is the reciprocal of the smallest k such that the system is unstable in 
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kJ$>. It is easy to see that the relation between the real SSV p and the real MSM kj^ is 

[i = 1 / k M ( 2 ‘ 48 ) 

As mentioned earlier, several significant results [17-22] have been obtained in the 
perturbed-characteristic-polynomials approach. Probably the most famous are the 
Kharitonov’s theorems [17] which apply to the special case with a hyper-rectangular 
perturbed region in the coefficient space. In this special case, the coefficients of the 
characteristic polynomial vary independently and the robust stability of the system can be 
easily determined by four bounding characteristic polynomials. Unfortunately, the 
Kharitonov's theorems cannot be applied to our problem since the coefficient variations of 
the characteristic polynomial are not independent. 

Bartlett, Hollot, and Lin [18] developed an important theorem which is applicable to 
the case when the coefficients of characteristic polynomial are linearly dependent. The 
theorem is now well known as the Edge Theorem: For the set of characteristic polynomials 
inside a polytope V in the coefficient space, every polynomial in (P is stable if and only if 
all the exposed edges of V are stable. This simplifies the stability checking tremendously 
since checking the stability of exposed edges is much simpler than checking that of the full 
IP. The exposed-edge stability checking is done by sweeping t from 0 to 1 such that 

ta‘+ (l-t)a j (2-49) 

are all stable for all vertices a 1 and of IP . 

Bialas [20] and Fu and Barmish [19] reduced the checking of the exposed-edge 
sweep stability to a one-shot test. They showed that t ft* + (1-t) cf is stable for all t e 
[0,1] if and only if the real eigenvalues of -HIT 1 are all negative where a 1 and a J are 
assumed to be stable and R and H. are the Hurwitz matrices for a 1 and a J respectively. 

Recently, a fast algorithm based on Chapellat and Bhattacharyya's Segment Lemma [22] 
was proposed by Bouguerra, Chang, Yeh, and Banda [23] for checking the stability of the 
exposed edges. The computation in the algorithm mainly depends on the number of vertices 
instead of the edges and therefore reduces the computation burden due to the "combinatoric 
explosion". 

There are no such celebrated properties in the parameter space as those in the 
coefficient space discovered by Kharitonov [17], Bartlett, Hollot, and Lin [18]. The 
closed-loop system may be unstable inside JS) although it is stable at all the edges of the 
hypercube <©. So far, there is no easy way of checking robust stability in the parameter 
space. 
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For each parameter vector 8 in the parameter perturbation domain J 9 , there is a 
corresponding characteristic polynomial, i.e., a coefficient vector a in the coefficient 
space. Let d(JS>) be the image of J9 in the coefficient space. The closed-loop system is 
robustly stable in J9 if and only if every characteristic polynomial in d ( JCD ) is stable. 
Although several significant results for robust stability have been obtained in the coefficient 
space, there is no efficient way to check robust stability for d(J9) since d ( J9) usually is 
neither a Kharitonov’s hyper-rectangle [17] nor a polytope considered by Bartlett, Hollot, 
and Lin [18]. 

Define the polytope iP(J§) as the convex hull of the 2 m image points in the (n+1)- 
dimensional coefficient space mapped from the vertices of JS) where n is the degree of the 
characteristic polynomial. If the mapping is multilinear then the image of the edges of the 
hypercube J9 will be the straight line segments connecting the corresponding mapped 
vertices. The image of <0, d(JD), is a subset of d(j§) and therefore is a subset of the 
polytope 1P(J9). Under the condition of the multilinear mapping, the stability of the edges 
of iP(J0) will guarantee the robust stability in J§. The multilinear mapping can be easily 
achieved by assuming that the nominal system M(s) in Fig.2.9 is a rational matrix with real 
coefficients. 

Now, we have an easy way to check the sufficient condition for the robust stability 
in JS) by using its corresponding polytope IP(JSj). The sufficient condition is still not 
enough to determine the real MSM k M . Any k such that the polytope IP(kJS)) remains 
stable, say k L , can be served as a lower bound for k M . However, there may exist some k > 

k L such that kJS> is stable although the corresponding polytope P(kJ§) is unstable. 

Any k which cause instability of k JS) or d(kJ9) qualifies as an upper bound for k M . 
To facilitate the description of the relations between the parameter space and the coefficient 
space, let the edges (vertices, resp.) of 1P(kJ$i) which are mapped from the edges (vertices, 
resp.) that are parallel to an axis of coordinates of kJ9 be called the crucial edges (vertices, 
resp.) and those which are not crucial be called noncrucial edges (vertices, resp.) . The 
noncrucial edges include two kinds of edges: su pplemental edges and fictitious edges . The 
supplemental edges are the image of the edges of kJS which are not in kJS). The fictitious 
edges are the edges of !P(kJ§) which are not mapped from the edges of kJS). The crucial 
edges are all in d(kJ!9). Thus, some k which causes instability at the crucial edges of 
P(kJS)), say k^, may be used as an upper bound for k M> If the lower and the upper bounds 
coincide or are close enough, we have the real MSM k M and the real SSV p. The objective 
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of the iterative algorithm to be presented in the paper is to reduce the gap of the lower and 
the upper bounds. When the gap is smaller than the desired accuracy e, i.e., Iky - kj < e, 

we have the real MSM k= k, and the real SSV u = 1/k, „. 

ML ^ M 

2.6. CONTROLLABILITY AND DEGREE OF CONTROLLABILITY 

The work performed in controllability has concentrated on discrete-time model 
representation of linear systems, specifically the deterministic discrete-time system model 
(sampled data system model) given by: 

i(k+l) = <t>(k) i(k) + 0(k) u(k) ; k = Icq, - (2-50) 

where x(k) is the system model state vector representation in R n at time k 
u(k) is the system model control vector in R m at time k 

<I>(k) is the (n-by-n) transition matrix of the system model from the state x(k) to 
x(k+l) with zero input applied to the system 
0(k) is the (n-by-m) input matrix of the system model at time k 
The degree of controllability measures the relation between the input u(k) and the state x(k) 
[28]. Recall that the system is not completely controllable if there exists some direction in 
state space which cannot be influenced by a control excitation. Further-more, some 
directions are more easily influenced than others. Thus, the Degree of Controllability 
(DOC) of a system must be some measure of the input-to-state gain along the direction in 
which it is smallest, that is the smallest gain of the input-state map. This gain vanishes if 
the system is uncontrollable. To define the DOC mathematically , consider the relevant 
state and input sets for the deterministic linear discrete-time system of Eq.(2-50): 

X := { x(k)€ R n , for ^>0, ke [k^ 1^+1, 1^+2 , )} = set of the states 

U := { u(k)€ R m , for ^>0, k€ [k^,, ko+1, kg +2 ) } = set of the inputs 

X and U are normed linear spaces with norm II • ll p of dimension n and m, respectively, 
where x(k)e X and u(k)e U. In terms of these sets, the DOC is defined as: 

Definition D.l (DOC): The degree of controllability, DOC(k N ,k 0 ) of a linear 

discrete-time system of Eq.(2-50) over N steps in the interval [k^.k^], is the minimum gain 
of the linear operator, G, at the initial state, Xo = x(ko). which describes the input-state 
relation of (2-50), that is [28] : 
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(2-51) 


f II Ax(k) II ) 

where G is the relation between the input space U and the state space X, namely: 

G = { (Ax(k),u(k)) : Ax(k)e X , u(k)e U , (Ax(k),u(k))*(0,0) & ( 00 , 00 ) } (2-52) 

Ax(k) = x(k) - x(k) is the solution to Eq.(2-50). 

Note that if the norm p=2 then the input space U is IMI 2 -space and the DOC can be 
considered as representing the degree of control accomplished when using an energy 
bounded controller. Moreover, if the system is not controllable then DOC = 0. 

To illustrate the mapping G of (2-52), we note that at k=k N , consider the solution of Eq. 
(2-50) given by: 


k - 1 

*(k) = <D(k, ko) ao + X F c (k- 1 . i) u(i) (2 -53) 

i-k Q 

where F c (k,i)= O(k,i)0(i), and <I>(k,i) is the transition matrix of the system from the state 
at time i to the state at k with zero input. This solution can be rewritten as: 

i(k N ) = F(k,ko) x(ko) + M c (k N ,ko) LL(k N ,ko) (2-54) 

where M c (k,ko) is an (n-by-mxk) controllability matrix given by: 

M c (k,k 0 ) = [O(k-l) | <J>(k,k-l)0(k-2) | d>(k,k-3) 0(k-3) | - |0(k,l) 0(k o )] (2-55) 

U(kjto) is a vector in R mx k made up of a sequence of input vectors given by: 

U(k,ko) = [ u T (k-l) I u T (k-2) I u T (k-3) I - I u T (ko) ] T (2-56) 

Thus the mapping G is the matrix M c (k N ,ko) of (2-55). This mapping is not 1:1, since Eq. 
(2-54) represents n equations in rrok unknowns, U(k N Jc()). This means that there are 
many control sequences (u(0), u(l), u(2),...,u(k-l)}that will return the system states to the 
origin at k N . The usual approach to solving Eq(2-54) for U(k N J cq) is to obtain the mean 
square solution (i.e, the pseudo inverse of M c (k N ,k 0 )), namely 

iI(k N ,ko) = - M c T (k N ,ko) ^(kovkjsj) [<I>(ko,k N ) M^k^ko) M c T (k N ,k 0 ) <J> T (k(),k N ) ]' 1 x(k 0 ) 
Which can be written using the controllability grammian W c (k N , ko) as: 

H(k N ,ko) = - M c T(k N ,ko) <D T (ko,k N ) [W c (k N , Icq) ]-l x(kq) (2-57) 

where W c (k N ,ko) is the (n-by-n) controllability grammian matrix given by: 
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W c (k N ,k 0 ) =F(ko, k N ) [M c (k N ,ko) M c T(k N , ko)] FT(ko, k N ) 


(2-58) 


or 

kf'M 

W c (k N , ko) = <fc(ko, i) 0 (i) 0 T (i)cD T (k o , i) ( 2 . 59 ) 

i = k 0 

Note that if the matrix W c (k N ,ko) is singular then the system is uncontrollable. The only 
way W c (k N ,ko) can be singular is if the matrix M c (k N ,ko) is not of maximal rank (i.e., 
rank{ M c (k N ,ko) }< n). This can be shown by observing that each column vector of 
M c (k N ,ko) represents a vector in the state space along which control is possible. For the 
system to be controllable in the whole space R n , it must be controllable in n linearly 
independent directions in the state space, or in terms of Eq.(2-54), the matrix M c (k N ,ko) 
must contain n linearly independent column vectors which form a basis in R n (that is, the 
range space of M c (k N ,ko) is equivalent to the state space, R n ). 

Next consider the evaluation of the DOC of Def. Dl. If the x s P ecif,ed is the origin (i.e., 
^specified = 0 ) then the Dx(k) = x(k) of Eq.(2-54) and the DOC of Eq.(2-51) can be 
evaluated using Eq.(2-54). Specifically, if the norm in Eq.(2-51) is p=2, then the N step 
degree of controllability at kq for the linear system of Eq.(2-50) is the square root of the 
smallest eigenvalue of the controllability grammian W c (k N , kq) [28], that is: 

DOC!(k N , kq) = [ X min (W c (k N , kq) )] 1/2 (2-60) 

where A 1tl , n ( W c (k N , kq)) is the smallest eigenvalue of the controllability grammian matrix 

Fact Fl(Muller and Weber [29] ): The degree of controllability of Eq. (2-60) represents 
the gain G of linear system of Eq.(2-50) under minimum energy , i.e., when 
Ilu(k)ll 2 = minimum , where 

m 

llli(k)ll 2 = [ X lui(k )|2 ] 1/2 
i=l 

Moreover, the direction in the X-space associated with A. m j n (dictated by the 
eigenvector of W c (k N , kq)) is the worst controllable directions of the system of 
Eq.(2-50), namely the influence of the control action is smallest as compare to the 
other directions. 

Remark: The DOC of Eq.(2-60) and Fact FI give only part of the information provided 
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by the controllability grammian , W c (k N , kg), about the controllability of system 
of Eq.(2-50). For example, if W c (k N , kg) is orthogonal, then according to Fact 
FI all direction in X-space are "worst controllable", so that one needs other 
measures which give structural information about the mapping G in addition to 
the DOC of (2-60). 

A viable candidate which gives additional structural information about the mapping G is the 
condition number of W c (k N , kg) [30]. This measure gives another DOC: 

Definition D2 (DOC based on Condition Number of W c (k N , kg) [30]): 

The degree of controllability of a deterministic linear discrete-time system, 
Eq.(2-50), is given by the inverse of the condition number of the 
controllability grammian matrix: 

1 11 W"i,(k N dc 0 ) h 

cond(W cN (k N ,k 0 )) IIW cN (k N ,ko)ll 2 (2 ' 61) 

where cond(A) is the condition number of matrix A 

iiaii 2 

cond(A) = — (2-62) 

IIA 

(Note that II A ll 2 is the largest eigenvalue, A. max , of A and II A - l \\ 2 is the 
smallest eigenvalue, A^, of A.) 

Fact F2: DOC 2 (k N ,k 0 ) takes on values in 0 < DOC 2 < 1, where DOC 2 = 0 implies that 
the system is uncontrollable while DOC 2 = 1 implies that all direction in the X- 
space are equally controllable in terms of requiring the same control energy (since 
the controllability grammian W c (k N ,kg) matrix is orthogonal). 

Additional structural information about the mapping G can be obtained by considering the 
Recovery Region, R(k N kg) [31-36], which is the set all the initial states, x(kg), that 
can be returned to the zero state in N steps by any bounded control, that is: 

Definition D3 (Recovery Region in terms of finite-time bounded-control): 

The N-step Recovery Region, R(N) at kg of a discrete-time system is 
defined as the set of initial states, xg, that can be returned to the origin in N 
steps in the interval [kg, k N ] with a bounded control: 
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R(N) = {xo:= x(ko),e R n I there exists a u(k)e U, ke [kQ,k N ], such that x(k N )=0 in 
N steps } 

where U := { u(k)e R m , luj(k) I < 1 for i=l,2,...,m }. (2-63) 

For the linear discrete-time system of Eq.(2-50), the recovery region, R(N), can be defined 
as (using the solution of Eq.(2-54), when x(k N )=0): 

R(N) = { xo e R n | xo = -[4>(k,ko)]-l M c (k N ,ko) U(k N ,ko), l Ui (k) I < 1 V i } (2-64) 

This R(N) is a polytope in the n-dimensional X-space as shown in Figure A- 1 for a 2nd 
order system. Note that the equation describing the R(N)-region is characterized by many 
constraints and only the binding ones are the relevant ones. 

Since in general, it is impractical to monitor the DOC using the R(N)-region 
characterization given in Eq.(2-64), we shall, instead, characterize it in terms "size", 
"shape" and "location", specifically: 

• The "size" is indicated by the volume of R(N)-region, Vol{R(N)}. 

• The "location" is indicated by the "center" of the R(N)-region, denoted by 

• The "shape" is indicated by the positive definite matrix Q which defines the largest 
imbedded ellipsoid. 

All the above structural information about the R(N)-region can be obtained by solving the 
following optimization problem: 

Problem P.l (Approximation of the Recovery Region, R(N»: 

Given the recovery region, R(N), of Eq-(2-64), find the largest approximating hyper- 


ellipse, E, contained in R(N), that is: 


DOC 3 (N) = Max { VOL(E)} 
E 

(2-65) 

where 


E = {*(0):ll [CKio-^lb^r) 

(2-66) 

B r 2 

V0L ®=,CI 

(2-67) 
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R i s the recovery region, R(N) 

C is an upper triangular matrix, whose diagonal elements Cjj > 0, for all i 
x c is the center of the recovery region 
r is the radius of the hyperellipse 
T(.) is the gamma function 

Remark: (i). The DOC3 of Eq.(2-65) represents the approximate"size" of R(N). 

(ii) The vector x c is the center of R(N). 

(iii) The matrix C gives the shape of R(N), since Q = C T C is the positive definite 
matrix defining the hyperellipse E. 

An illustration of the solution of Problem P.l is shown in Figure A.3 which show the 
largest hyperellipse embedded in the R(N)-region for a 2nd order system. 

Fact F3: The radius r of the hyperellipse Eof Eq. (2-66) gives an approximate value of 
the gain G of linear system of Eq.(2-50). Moreover, the eigenvectors associated 
with A. m i n of C-matrix gives the worst controllable directions of the system of 
Eq.(2-50), namely the direction in which the influence the control action is 
smallest as compare to the other directions. 

A set of software was developed to solve the optimization problem of Problem PI. It is 
described in Appendix A. This software investigates the robustness of linear time invariant 
system, and determines the degree of controllability of the system under uncertainty d. The 
function of the software is: 

1. To solve the optimization problem Problem PI. 

2. To graphically display the R(N)-region in X-space of the LTI system. 

3. To graphically display the regions of 6-space of the LTI system. 

4. To display the approximating hyperellipse to the R(N)-region in the X-space. 

5. To display the approximate regions in the 6-space which are based on a 
prescribed degree of controllability (DOC). 

Figure 2.10 gives a structured diagram of this software and shows some of the specific 
capabilities of the software. There are two main options of the controllability robustness 
software One, investigates the recovery region while the other investigates the allowable 
uncertainty region in the d-parameter space. Referring to Fig. 2.10, the branch denoted by 
SYSDAT is the data file containing all the information about the LTIS and the parameter 
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Fig.2.20 Software structure for studying the controllability robustness of LTI systems 






















uncertainty. The branch whose heading is the Exact Recovery Region obtains the 
R(T=NAt)- region of Definition D3, where T=NAt is the total time of the application of the 
control. Specifically, the software computes the set of linear inequalities { Fx < b } and 
determines the binding constraints which define R(NAt). This is illustrated in Fig. 2.1 1, 
which shows the region R(NAt) when T = NAt = 0.3 sec. for a 2nd order linear system 
given by: 


or 


G(s) 1 


s 2 + 3s + 2 


1 

o 

1 

1 

X 

+ 

1 

L-2 -3J 

L x 2 J 

l 


U 



Fig-2. 1 1 : The recovery region for the 2nd order system, T=0.3 


The discrete time system representation using Eq.(A-3) of the appendix when At = 0.075 
sec and n= 15 , is the system model of Eq.(2-50) when 


0.995 .067 
.-.134 .794 J ’ 



.0026 
.067 . 


Referring to Fig. 2.10, the branch whose heading is the Approximate Recovery 
Region, AR(T), finds the best approximation of the R(T=NAt)-region given by the 
hyperellipse E of Eq. (2-66). Specifically, it solves Problem PI to obtain the shape (matrix 
C), the location (center x c ) and the size (radius, r) of the R(T=NAt)-region. Figure 2.12 
illustrates the AR(T=NAt)-region for the R(T=NAt)-region shown in Fig. 2.11. The 
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radius r = 0.079 and the matrix C is 


C = [4-697.6834 
. 0 .2129. 



Fig.2. 12: The approximate recovery region for the 2nd order 
system when the control period T=0.3 

Referring to Fig. 2.10, the branch whose heading is the Allowable Uncertainty 
region AA(T), finds the set of values of the uncertainty parameters in the 8-space. 
Specifically, 


AA(T-NAt) - { 8 : G(s;8) gives a level of controllability as prescribed by a 
R(T=NAt)-region } (2-69) 

where G(s;8) is the transfer function of the LTI system. This AA(NAt) can be defined in 
terms of the structural indicators of R(NAt) given by the C-matrix, center x c and radius r as 
follows: 


AA(NAt) - { 8 : G(s;8) gives a level of controllability specified by triplet (C,x c ,r) } 

(2-70) 

At present, this part of the software has not yet been developed. What has been 
implemented are the two branches shown in Fig. 2.10 under the heading Recovery 
Region. Based on preliminary tests, the software works well except when the 
hyperellipse is ill-conditioned (i.e., when the C-matrix is ill-conditioned). We suspect that 
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this is due to our choices of the termination conditions in the (C,r) optimization algorithm 
of Fig. A. 3 and the overall optimization method of Fig. A.4 given in Appendix A. 

To initiate the design of the software indicated by branch in Fig. 2.10 whose heading is the 
Allowable Uncertainty region, AA(T), we performed sensitivity studies to determine 
the type of variation to be expected in R-region when the system parameters change. 
Specifically, for the system 


G(s) = — — 1 

s 2 +as + h 

we studied the effect of variation of the coefficients a and b on the triplet (C,x c ,r). Figure 
2.13 shows the R-region when a = b = 1.414. Note that this R-region changed both in 
size and orientation from that shown in Fig. 2.11. This changed is measured by a different 
values of r = .0464 and C-matrix 

c _ T 4.76 0.8 ' 

0 .21 . 

The center x c = 0, the origin of the X-space. This is the usual case because of the 
definition of the R-region. 



Fig.2.13: Recovery region for a 2nd order system, T=0.3 


38 




3. SUMMARY OF THE WORK 

Task 1: Design of Robust Optimal Controllers via H~ Approach 

Task 1.1: Study of the aircraft flight guidance/control models and problems 
provided by the Aircraft Guidance and Control Branch, NASA Langley 
Research Center. 

Status: In the process of obtaining appropriate practical aircraft flight 
guidance/control models from AGCB scientists. Before the models arrive, 
we use a hypothetical model with plant perturbations and uncertain 
disturbances. 

Accomplishment: A procedure to represent plant perturbations and uncertain 
disturbances as H°° models is presented. 

Task 1.2: Formulation of the aircraft flight guidance/control problems as a robust 
H“ optimization problem with unstructured and parametric uncertainties. 

Status: Before the practical models arrive, a hypothetical model is used to 
demonstrate how to formulate a robust H“ optimization problem. 

Accomplishment: A detailed procedure to formulate a robust H~ optimization 
problem is presented. See [Pub 1,2] for details. 

Task 1.3: Computation of the H°° norm of a given transfer function. 

Status: 100% of the proposed work has been done. See [Pub 3] for details. 

Accomplishment: The pole-zero diagram of G(s), the root locus concept, and 
the jco-axis transmission zeros of y 2 I-G T (-s)G(s) are employed to quickly 
locate a narrow frequency interval which brackets the supremum of 
o[G(j(«>)]. Then Brent's method (an improved parabolic interpolation 
search method) is employed to search for the supremum of o[G(jco)] over 
the frequency interval. 

Task 1.4: Solution of the robust H“ optimization problem. 

Status: 100% of the proposed work has been done. See [Pub 4,5,6] for details. 

Accomplishment: Some useful properties of the two H°° Riccati solutions have 
been discovered. Among them, the most prominent one is that the spectral 
radius of the product of these two Riccati solutions is a continuous, 
nonincreasing, convex function of y in the domain of interest. Based on 
these properties, quadratically convergent algorithms are developed to 
compute the optimal H°° norm. 

Task 1.5: Algorithm development for the computational problems arising in the 
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robust H°° optimization problem. 

Status: 60% of the proposed work has been done. See [Pub 4-8] for details. 

Accomplishment: The algorithm to compute the optimal H~ norm and the 
construction of an optimal or a suboptimal H~ controller are presented. 

Task 1.6: Computer simulation. 

Status: 40% the proposed work has been done. 

Accomplishment: The time-domain and frequency-domain behavior of the 
closed-loop system for the hypothetical model has been tested. 

Task 2: Design of Robust Optimal Controllers via Mixed H 2 /H°° Approach 

Task 2.1: Formulation of the aircraft flight guidance/control problems provided by 
the Aircraft Guidance and Control Branch, NASA Langley Research 
Center as a mixed H~/H optimization problem with unstructured 
uncertainties. 

Status: In the process of obtaining appropriate practical aircraft flight 
guidance/control models from AGCB scientists. Before the practical 
models arrive, a hypothetical model is used to demonstrate how to 
formulate an H^/H 2 optimization problem. 

Accomplishment: A procedure to formulate an H°°/H 2 optimization problem is 
presented. See [Pub 9-12] for details. 

Task 2.2: Solution of the mixed H 2 /H°° optimization problem. 

Status: An important special case of the mixed H 2 /H~ optimization problem has 
been completely solved. 

Accomplishment: In case that the H 2 and the H°° transfer functions are identical, 
the mixed H 2 /H°° problem can be solved by H°° technique with y a free 
parameter. Then y is used to trade off H~ robustness against H 2 
performance. See [Pub 9-12] for details. 

Task 2.3: Plotting the H 2 /H“ curve and using it for robust controller design. 

Status: 80% of this study has been done. 

Accomplishment: A graphical method to determine the best y for H" robustness 
and H 2 performance is proposed. See [Pub 1 1,12] for details. 

Task 2.4: Algorithm development for the computational problems arising in the he 

mixed H /H“ optimization problem. 

Status: 40% of the proposed work has been done. See [Pub 1 1,12] for details. 
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Accomplishment: Algorithms to generate the H 2 /H°° curve and the construction 
of an H 2 /H°° controller are presented. 

Task 2.5: Computer simulation. 

Status: 40% the proposed work has been done. 

Accomplishment: The time-domain and frequency-domain behavior of the 
closed-loop system for the hypothetical model has been tested. 

Task 3: Construction M-A Structure 

Task 3.1: Construction of an M-A structure of a system with parametric 
uncertainties for robust stability analysis. 

Status: 100% of the proposed work has been done. See [Pub 16-17] for details. 
Accomplishment: A systematic way of constructing an M-A structure of a 
system with parametric uncertainties for robust stability analysis is 
presented. 

Task 3.2: Construction of an M-A structure of a system with parametric and 
unmodelled uncertainties for robust stability analysis. 

Status: 100% of the proposed work has been done. See [Pub 16-17] for details. 
Accomplishment: A systematic way of constructing an M-A structure of a 
system with parametric and unmodelled uncertainties for robust stability 
analysis is presented. 

Task 3.3: Construction of an M-A structure of a system with parametric and 

unmodelled uncertainties for robust stability and robust performance 
analysis. 

Status: 100% of the proposed work has been done. See [Pub 16-17] for details. 
Accomplishment: A systematic way of constructing an M-A structure of a 

system with parametric and unmodelled uncertainties for robust stability 
and robust performance analysis is presented. 

Task 3.4: Verification of the minimality of an M-A structure of a system. 

Status: 30% of the proposed work has been done. See [Pub 16-17] for details. 
Accomplishment: We show that if the dimension of A is the least number of 

assigned parameters such that the coefficients of the corresponding 
characteristic polynomial are multilinear functions of the assigned 
parameters, then the M-A structure is minimal. 



Task 3.5: Construction of a minimal M-A structure. 

Status: 30% of the proposed work has been done. See [Pub 16-17] for details. 

Accomplishment: For the case that the state-space realization is a linear function 
of the uncertainties, a minimal M-A structure can be constructed in a 
systematic way. 

Task 4: Controllability and Observability of Perturbed Plant 

Task 4.1: Controllability. 

Status: 70% of the proposed work has been done. 

Accomplishment: Developed (i) a software package for investigating the 
controllability robustness of linear time invariant systems. This software 
obtains the size, shape and location of the recovery region, and (ii) 
relation between the various DOC measures encountered in the literature. 

Task 4.2 Observability 

Status: No work performed. 

Task 4.3 Integration of Controllability and Observability Concepts into 
the Overall Robust Controller design 

Status: No work has been performed as yet. 
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4. CONCLUSION AND WORK FOR FUTURE RESEARCH 
4.1 Conclusion 

In [Pub 4] we revealed several important eigen properties of the stabilizing 
solutions of the two H“ Riccati equations and their product. Among them, the most 
prominent ones are: (1) plX^ylY^y)] is a nonincreasing, convex function of y on ((3,+o°). 

(2) A. min [&(y)] is a nondecreasing, concave function of yon (ct^+oo). (3) &(y) is invertible 
almost everywhere on (a x ,+°°). Based on these properties, quadratically convergent 
algorithms [Pub 5,6] were developed to compute the y^, the optimal H°° norm, such that 

pfXeofyJY^fy.,)] = yf and to compute (3, the infimum of y such that the two H°° Riccati 
equations have positive semidefinite stabilizing solutions. 

The formulation of the standard H~ optimization problem, an easy way of 
constructing a state-space realization of the generalized plant, an efficient algorithm for 
computing the optimal H“ norm, and a modified version of Glover and Doyle formulas for 
constructing an optimal controller were addressed in the paper [Pub 2]. No numerical 
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difficulty will arise in constructing an H°° optimal controller if we are allowed a proper 
controller. In most applications, we may like to have a strictly proper controller with limited 
bandwidth. In this case, a trade-off between the H“ performance and the bandwidth should 
be made by degrading the H°° norm from its optimum in order to reduce the controller 
bandwidth. 

In [Pub 11], we found an easy way to trade off the H°° robustness and H 2 
performance for the case that H°° and H 2 transfer functions are the same. In general, these 
two transfer functions may be different. In [Pub 12], we formulated a more realistic control 
problem with H 2 performance index and H°° stability robustness constraint into a mixed 
H 2 /H°° problem. No optimal solution yet is available for this more general mixed H 2 /H°° 
problem. Although the optimal solution for this mixed H 2 /H°° control has not yet been 
found, we proposed a design approach which can be used through proper choice of the 
available design parameters to influence both robustness and performance. 

The theoretical work in controllability on developing a relation between the various 
degree of controllability (DOC) measures encountered in the literature, indicates that there is 
a relation and what is needed now to establish the specific relations. The simulation work 
concentrated on developing software for investigating the controllability robustness of 
linear time invariant systems. A preliminary software package (described in Appendix A) 
was developed for measuring the size, shape and location of the recovery region. Initial 
indications show that the various optimization algorithms (needed to automatically obtain 
the value of the indicators) work, however they may need some refinements in order to take 
care of such problems as obtaining DOCs for ill-conditioned recovery regions. 

4.2 Work for Further Research 

4.2.1 M-A Structure and Robust Stability Analysis for Structured 
Uncertainties 

M-A structure is a rearrangement of a perturbed system where M(s) is the nominal 
system and A is a block diagonal matrix which consists of all the perturbations. M-A 
structure is essential in the SSV(structured singular value) analysis and design techniques. 
Although A is assumed stable, the M-A structure can handle a large class of uncertainties 
which include unstable perturbation, i.e., the number of unstable poles of the perturbed 
plant can be different from that of the nominal plant. It is always possible to pull out all 
uncertain parts from a perturbed system and form an M-A structure. However, very little 
about how to obtain this structure has been addressed in the literature. A systematic 
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procedure for constructing M-A structure is proposed. 

A minimal M-A structure means that the dimension of A (or M) of the M-A 
structure is minimal. We can construct an M-A structure for a given perturbed system, but 
the dimension of the structure may be unnecessarily large. A nonminimal structure will 
cause unnecessary complexity in computation and therefore a minimal M-A structure is 
essential in the computation of the SSV. So far, only some special cases of the minimal M 
A structure problems have been addressed [37,38]. A progress in the construction of a 
minimal M-A structure will greatly simplify robust stability analysis for structured 

uncertainties. 

4.2.2 Design of Robust Flight Control Systems via H~ Optimization 

In this research period, we developed a fast algorithm to compute the optimal H~ 
norm and to construct an optimal H~ controller without numerical difficulty. Next, we will 
apply this design tool to flight control systems. We will use the H norm of the 
complementary sensitivity function as a measure of robust stability and will formulate an 
H°° optimization problem which minimizes the maximum error energy subject to a robust 
stability constraint. The error reduction and the robust stability can be traded off by 
choosing the weighting matrices in the cost function. Initially, a weighting matrix is chosen 
and the corresponding H~ optimization problem is solved to obtain an optimal controller. 
The weighting matrices are modified iteratively until the robust stability constraint is just 
satisfied. The way the weighting matrices affect the trade-off is only partially understood. 
In addition to the magnitude of the weighting matrix, the structure of the weighting matrix 
is also an important factor in the trade-off. We will investigate how the weighting matrices 
affect the trade-off. We will also develop an iterative updating procedure for the weighting 
matrices by which a robust controller can be designed such that the closed-loop system is 
robustly stable and the maximum of error energy is minimized. 

4.2.3 Design of Robust Controllers via Mixed H 2 /H~ Optimization 

The optimal H~ controller which minimizes the H~ optimization problem usually 
has wide bandwidth and leads to a poor H- performance. In [13], we found that a little bit 
sacrifice of the H~ norm will greatly reduce the bandwidth of the controller and improve the 
H 2 performance tremendously. It is practical to formulate a robust control problem as a that 
of minimizing an H 2 cost function subject to an H~ bound. The H 2 cost function is 
equivalent to the well known LQG cost function and the H°° bound can take care of robust 
stability and robust error reduction. 
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In [39], we formulated a more realistic control problem with H 2 performance index 
and H°° stability robustness constraint into a mixed H 2 /H°° problem. The H 2 and H°° 
transfer functions are different. No optimal solution yet is available for this more general 
mixed H 2 /H°° problem. A research in finding an efficient approach for the mixed H 2 /H“ 
optimization problem is proposed. 

4.2.4 Controllability and Observability of Perturbed Plant 

For task 4 in the next research period we will concentrate on: (i) Theoretical work 
on developing a relation between the various DOC measures encountered in the literature, 
and determine which of these is most suitable for characterizing the Allowable 
Uncertainty region, AA(T), discussed in Section 2.6. (ii) Improving the algorithms 

used in developing the Recover Region software shown in Fig. 2. 10 Specifically, we will 
take care of the ill-conditioning problem currently encountered when the recovery region is 
very narrow. We will also design and develop the Allowable Uncertainty region, 
AA(T), software indicated in Fig. 2.10. (iii) Perform robustness studies for a variation in 
the system parameters (c,} and [bj } of the system function 


G(s) 


XbfU 
Z qsi 
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Appendix: Controllability Robustness Software 

A.l Purpose 

The purpose of the Controllability Robustness (COR) software is to investigate the 
effectiveness of the various performance indices used for measuring the degree of 
controllability of linear time invariant system under uncertainty 5. The function of the 
software is: 

1. To graphically display the Recovery region, R(T), in the X-space. 

2. To automatically compute the degree of controllability (DOC) indices which are needed 
to characterize the controllability of the LTI system. 

3. To graphically display the Approximate Recovery Region, AR(T) in the X-space. 

4. To graphically display the set of parameter variations in the 8-space which are 
guaranteed to give a specified level of system controllability. 

The graphical capabilities is illustrated in Fig. 2.10 (of the report ) which gives a structured 
diagram of this software. 

A.2 Method 

The method used in obtaining the DOC indices of system is summarized in the flowchart of 
Figure A.l. The inputs to the program are the matrices {A,B} and the parameter 
variations, d, of the linear time invariant system 

x = A(8)x + B(5)u (A . 1} 

Remark: It is planned that once a method for the M-A decomposition will be available, 
then the input to the COR program will be the G(s,8) system function. 

The system of Eq.(A-l) is sampled with a sampling period. At, to obtain the linear discrete 
time model of Eq.(2-50), namely 

x(k+ 1 ) = 0(8) x(k) + 0(8) u(k) (A-2) 

where the matrices O and 0 are obtained using the power series expansion of the 
transition matrix 

O = [ I + AAt + (AAt) 2 /2! + ... + (AAt) n /n! ] (A-3a) 

0 = [ I + AAt + (AAt) 2 /2! + ... + (AAt) n /n! ] At (A-3b) 

This is done in the subroutine DSCRET shown in Fig. A.l. Referring to Fig. 2.10, the 


50 



program which computes the Recovery Region, IR(T=NAt), implements the definition 
of Eq.(2-64) and computes a set of linear inequalities { Fx < b} and then identifies those 
constraints which are the binding constraints of R(NAt). This is implemented in subroutine 
CONTRG shown in Fig. A.l. 

The Approximate Recovery Region, AR(T), software solves the optimization 
Problem PI to obtain the shape (matrix C), the location (center x c ) and the size (radius, r) 
of the R(T=NAt)-region and then uses this to define the largest embedded hyperellipse in 
the R(T=NAt)-region. The subroutine which executes the solution of Problem PI is 
labeled ELLIPSE in Fig. A.l. The approach to solving Problem PI is based on the 
realization that when the matrix C is given then Problem PI is a Linear Problem (LP) 
problem. This is exploited in our method as summarized in Fig. A.2. The first step is to 
set up the linear constraints. Then, depending on if we want to obtain a new matrix C or 
not we solve either an Linear Problem (LP) problem or a Non Linear Problem (NLP), that 
is: 

LP problem finds the optimum center, x c , and radius, r, for a given matrix C. The 
LP problem solves a volume maximization problem in which the embedding hyperellipse E 
shape is held constant but its size (i.e., volume of Eq.(2-66)) is maximized by increasing 
the radius, r, and moving it around inside R (by moving the center, x c ) while satisfying the 

constraint ECR. 

NLP problem finds the matrix C and radius r for a given center x c . It solves the 
same volume maximization problem of Eq.(2-66), when the center x c is fixed. The C- 
matrix and r are optimized by a sequential unconstrained minimization (SUM) method. The 
method minimizes an objective function 7k which is a pill norm of a vector representing the 

slack variable of the inequalities. The details of the updating method is shown in Fig. A. 3. 
The optimum is found by computing the gradient of 7k with respect to tjj, the components 

of the inverse of the C matrix (T=C'l). (The reason for using T instead of C in the 
optimization is because the SUM method has better convergent properties.) 

The method for optimizing all three sets of parameters, i.e. C-matrix, x c and r is 
shown in Figure A. 4. Note that it uses both the NL and the NLP algorithms to solve the 
optimization problem stated in Problem Pl. The subroutine implementing the optimization 
algorithm shown in Fig. A.4 is called ELLIPSE. 




Figure A-l : Flowchart of the COR program for studying the Controllability Robustness of Linear 
Time Invariant Systems 
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Figure A. 2: Method for finding the largest hyperellipse, E, embeded in 
a polytope R(T) 
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T-MATRIX OPTIMIZATION 



Figure A. 3: Flowchart of algorithm for optimizing C- Matrix 
(Note - T matrix is the inverse of the C matrix) 
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Fig. A.4: Parameter optimization algorithm for computing C-Matrix, center x and radius r 
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